Reading

CIFTI files organize the gray matter of the brain into “grayordinates”: vertices representing the left and right cortical surfaces, and voxels representing the subcortical gray matter structures and the cerebellum. A CIFTI file consists of two parts: (1) an XML header which contains all the metadata including medial wall locations, subcortical structure labels, and the subcortical volumetric mask; and (2) a NIFTI format matrix representing all the grayordinate data. These components are read in together with read_cifti:

basename(cifti_fnames["dtseries"])

## [1] "Conte69.MyelinAndCorrThickness.32k_fs_LR.dtseries.nii"

xii <- read_xifti(cifti_fnames["dtseries"])

The result of read_cifti is a "xifti" object. We can get an overview of it using its summary S3 method:

xii # same as `summary(xii)`

## ====CIFTI METADATA===================
## Intent:           3002 (dtseries)
## - time step       1 (seconds)
## - time start      0
## Measurements:     2 columns
## 
## ====BRAIN STRUCTURES=================
## - left cortex     30424 data vertices
##                   2068 medial wall vertices (32492 total)
## 
## - right cortex    30527 data vertices
##                   1965 medial wall vertices (32492 total)

By default, read_cifti will read-in all existing brainstructures. Set brainstructures=c("left","right"), for example, to read in just the left and right cortex data. The full set of choices for brainstructures is any combination of "left", "right" and "subcortical", "all" for all three, or "existing" (default).

"xifti" objects are lists with entries data (the grayordinate data matrix, separated by brainstructure), surf (surface geometry), and meta (metadata, most of which is from the NIFTI XML header). surf distinguishes a "xifti" from a CIFTI file: the left and right cortical surface geometries are not included in CIFTI files, so they must be read from separate surface GIFTI files (ending in surf.gii). The surface must be compatible: the number of vertices must be the same, and each vertex in the CIFTI data must correspond to the vertex location in the corresponding GIFTI surface file. In this way, a "xifti" represents a combination of a CIFTI file with compatible GIFTI files for the cortical mesh.

We can add GIFTI surface geometries with add_surf:

xii <- add_surf(xii, surfL=surfL_fname, surfR=surfR_fname)
xii

## ====CIFTI METADATA===================
## Intent:           3002 (dtseries)
## - time step       1 (seconds)
## - time start      0
## Measurements:     2 columns
## 
## ====BRAIN STRUCTURES=================
## - left cortex     30424 data vertices
##                   2068 medial wall vertices (32492 total)
##                   left surface geometry is present
## 
## - right cortex    30527 data vertices
##                   1965 medial wall vertices (32492 total)
##                   right surface geometry is present

Alternatively, we could have provided the surface geometries at the outset of reading the CIFTI file:

xii2 <- read_xifti(cifti_fnames["dtseries"], surfL_fname=surfL_fname, surfR_fname=surfR_fname)
all.equal(xii, xii2) # same result

## [1] TRUE

To only read the CIFTI header, use info_cifti. Let’s read the header of the dscalar file:

xii_info <- ciftiTools::info_cifti(cifti_fnames["dscalar"])
str(xii_info, nchar.max=50) # shows header structure

## List of 3
##  $ cortex :List of 2
##   ..$ medial_wall_mask:List of 2
##   .. ..$ left : logi [1:5762] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. ..$ right: logi [1:5762] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ resamp_res      : NULL
##  $ subcort:List of 4
##   ..$ labels     : NULL
##   ..$ mask       : NULL
##   ..$ trans_mat  : NULL
##   ..$ trans_units: NULL
##  $ cifti  :List of 4
##   ..$ intent         : num 3006
##   ..$ brainstructures: chr [1:2] "left" "right"
##   ..$ names          : chr [1:2] "MyelinMap_BC_decurv" "corrThickness"
##   ..$ misc           :List of 4
##   .. ..$ ParentProvenance : chr "C:\\Users\\damon\\AppData\\Local\"| __truncated__
##   .. ..$ ProgramProvenance: chr "Connectome Workbench\nType: Comma"| __truncated__
##   .. ..$ Provenance       : chr "C:\\Users\\damon\\Desktop\\WORKBE"| __truncated__
##   .. ..$ WorkingDirectory : chr "C:/Users/damon/Desktop/ciftiTools/vignettes"

To read in only certain columns of a CIFTI file, use the idx argument:

read_xifti(cifti_fnames["dtseries"], idx=2) # second column only

## ====CIFTI METADATA===================
## Intent:           3002 (dtseries)
## - time step       1 (seconds)
## - time start      0
## Measurements:     1 column
## 
## ====BRAIN STRUCTURES=================
## - left cortex     30424 data vertices
##                   2068 medial wall vertices (32492 total)
## 
## - right cortex    30527 data vertices
##                   1965 medial wall vertices (32492 total)

Writing

When a "xifti" object is written to files, the CIFTI components are placed in a CIFTI file and the surface geometries, if any, are placed in GIFTI files.

out_dir <- "output"

write_xifti(
  xii, 
  file.path(out_dir, "my_cifti.dtseries.nii"), 
  file.path(out_dir, "my_L.surf.gii"), file.path(out_dir, "my_R.surf.gii")
)

## Writing left cortex.
## Writing right cortex.
## Creating CIFTI file from separated components.
## Writing surface geometry GIFTI(s).

With separate_cifti, a CIFTI can be separated and written into its component parts: the cortical data can be written to GIFTI metric or label files, and the subcortical data can be written to a NIFTI file. In addition, any ROIs or labels will also be written to files. The files are automatically named unless a new file name is provided.

# Use default names for everything except left cortex
separated_fnames = separate_cifti(
  cifti_fnames["dscalar_ones"], brainstructures="existing", 
  cortexL_fname="my_left_cortex.func.gii", write_dir = out_dir
)
# Files written to `out_dir`, or current working dir. if not specified
basename(separated_fnames)

## [1] "ones_1k.my_left_cortex.func.gii" "ones_1k.sep.ROI_L.func.gii"     
## [3] "ones_1k.sep.R.func.gii"          "ones_1k.sep.ROI_R.func.gii"     
## [5] "ones_1k.sep.nii"                 "ones_1k.sep.labels.nii"         
## [7] "ones_1k.sep.ROI.nii"

Separated files can be read into R with the oro.nifti, RNifti, and gifti packages, and combined into a "xifti" object with as.xifti.

Surface visualization

view_xifti_surface(xii) displays the cortical data on the surface mesh. This function has several primary arguments:

• color_mode specifies the nature of the data values: "sequential", "qualitative" and "diverging". If it is not provided, a default mode that makes sense for the data will be used.
• colors specifies the color palette to use. If it is not provided, a default palette that makes sense for the color_mode is used.
• idx controls which column(s) to display.
• widget and fname control the output type. If fname is not provided, an interactive plot is created: by default, an OpenGL window if the length of idx is one, and an embedded HTML widget if the length of idx is greater than one. widget can be used to override this default. On the other hand, if fname is provided, static image files (png) for each idx are created, unless fname ends in .html in which case an interactive html file will be saved. Lastly, both OpenGL windows and HTML widgets can be embedded in R Markdown documents for knitting; refer to the source code of this vignette to see how this works.
• surfL and surfR specify the surface geometry to plot the data on. If not provided, the surfaces in the "xifti" object is used. But if those are also unavailable, the “inflated” surfaces included in ciftiTools are used.

Let’s see an example using each color_mode option. Note how the included surfaces are used in the first plot, but if none are present as in the second and third plots, the default surfaces are automatically used for visualization. (In our case, the included surfaces are the default surfaces.) We’ll also make the second plot interactive by requesting display of two idx. Try clicking and dragging around the second plot to rotate, and scrolling to zoom in and out. Note that the first and third plots are OpenGL window snapshots, and that the second plot is an embedded HTML widget.

# Normally `cex.title` doesn't need to be set, as it defaults to a good choice.
#   But when knitting static images this way, the default becomes a bit too big
#   based on how knitting works.
view_xifti_surface(xii, idx=1, zlim=c(1,2), title='color_mode = "sequential"', cex.title=1.3)

dtseries file; first column; sequential palette

xii <- read_cifti(cifti_fnames["dscalar"]) # no GIFTI included, so the default inflated surface is used.
view_xifti_surface(
  xii, idx=1:2, zlim=c(0,5), color_mode = "diverging",
  title='color_mode = "diverging"', cex.title=1.3
)

view_xifti_surface(
  read_cifti(cifti_fnames["dlabel"]), 
  # Interactively, a color legend that displays the label names will also be printed.
  legend_ncol=5, 
  title='color_mode = "qualitative"', cex.title=1.3
)

dlabel file; first label; palette from label metadata

Volume visualization

view_xifti_volume(xii) displays the subcortical data in slices. To view interactively in an interactive session, set interactive=TRUE. By default, a series of slices is displayed overlaid on the MNI template. The orientation and numbers of slices can be adjusted. A .png or .pdf file can be written with fname. The same color arguments from view_xifti_surface work too: colors, color_mode and zlim.

# cifti_fnames["dscalar_ones"] is the only file with subcortical data
xii <- read_cifti(cifti_fnames["dscalar_ones"], brainstructures="subcortical")
view_xifti_volume(xii)

## `zlim` not provided: using color range 0 - 1 (data limits: 1 - 1).

## Selected axial slices: 9, 13, 17, 21, 25, 29, 33, 37, 41

Subcortical data (all ones)

# For information only, since papaya viewer cannot be embedded in knitted file
view_xifti_volume(xii, interactive = TRUE)

The S3 method plot(xii) will display all the data present in the xifti using view_xifti_surface, view_xifti_volume, or both. If both plots are made, the color mode and range will be shared between the two plots.

Tips and tricks

NA values will be uncolored. You can use this to, for example, only color values meeting a certain threshold.

xii <- read_cifti(cifti_fnames["dscalar"])
# Convert to z-scores
xii <- scale_xifti(xii)
# Threshold
xii <- transform_xifti(xii, function(x){ifelse(x<2, NA, x)})
view_xifti_surface(
  xii, title='MyelinMap_BC_decurv: z > 2', 
  cex.title=1.3, zlim=c(2,3), NA_color="#505560"
)

dscalar; using NA values to selectively color

Something to note is that the shadows from 3D rendering will darken colors in some folded regions of the brain. Another consideration is that the colors are smoothed continuously. In some situations, these default 3D shading properties may be undesirable. The shadows and material arguments can be used to adjust aspects of 3D shading.

view_xifti_surface(
  xii, title='MyelinMap_BC_decurv: z > 2, less shadows', 
  cex.title=1.3, zlim=c(2,3), NA_color="#505560",
  shadows=0
)

Lessening the shadows

view_xifti_surface(
  xii, title='MyelinMap_BC_decurv: z > 2, exact colors', 
  cex.title=1.3, zlim=c(2,3), NA_color="#505560",
  material=list(lit=FALSE, smooth=FALSE)
)

Exact coloring, by disabling shadows and color smoothing

Resampling

ciftiTools can resample CIFTI files to a lower resolution. Here, we resample the 32k dtseries file to 2k. (The number refers to the count of vertices on a single hemisphere.) We also provide the surfaces and resample them in conjunction.

resampled_xii_fname <- "my_new_resampled.dtseries.nii"
resampled_surfL_fname <- "my_resampled_surfL.surf.gii"
resampled_surfR_fname <- "my_resampled_surfR.surf.gii"
  
xii_2k <- resample_cifti(
  cifti_fnames["dtseries"], resampled_xii_fname,
  resamp_res = 2000,
  surfL_fname, surfR_fname,
  resampled_surfL_fname, resampled_surfR_fname,
  write_dir = out_dir
)

## Separating CIFTI file.
## Time difference of 0.3485959 secs
## Resampling CIFTI file.
## Time difference of 2.359837 secs
## Merging components into a CIFTI file... 
## Time difference of 0.03574109 secs

basename(xii_2k)

## [1] "my_new_resampled.dtseries.nii" "my_resampled_surfL.surf.gii"  
## [3] "my_resampled_surfR.surf.gii"

Resampling can also be performed while reading a CIFTI file into R.

read_cifti(cifti_fnames["dscalar"], resamp_res=2000)

## ====CIFTI METADATA===================
## Intent:           3006 (dscalar)
## - names           "MyelinMap_BC_decurv", "corrThickness"
## Measurements:     2 columns
## 
## ====BRAIN STRUCTURES=================
## - left cortex     1856 data vertices
##                   106 medial wall vertices (1962 total)
## 
## - right cortex    1863 data vertices
##                   99 medial wall vertices (1962 total)

Smoothing

Use smooth_cifti to perform smoothing. Like resampling, this function works on both CIFTI files and "xifti" objects.

smoothed_xii_fname <- "my_smoothed_cifti.dtseries.nii"

# Smoothing a CIFTI file
smooth_cifti(
  cifti_fnames["dtseries"], file.path(out_dir, smoothed_xii_fname),
  surf_FWHM=2, vol_FWHM=2,
  surfL_fname=surfL_fname, surfR_fname=surfR_fname,
  subcortical_zeroes_as_NA=TRUE
)

# Visualizing the smoothed file.
# Let's demonstrate the ability to use RColorBrewer palettes!
plot(
  read_cifti(file.path(out_dir, smoothed_xii_fname)), 
  surfL=surfL_fname, surfR=surfR_fname, 
  zlim=c(1,2), color_mode="diverging", colors="Spectral"
)

Smoothed CIFTI

Mesh operations

ciftiTools includes several functions for operations that depend on the vertex adjacencies of the surface mesh.

edit_mask_surf can erode, dilate, or extract the border vertices of a binary mask. (Here we’ll also demonstrate how to use view_comp to display multiple plots simultaneously, without the interactivity of a widget.)

xii <- load_parc("Schaefer_400")
xii <- remove_xifti(xii, "cortex_right")
parc_name <- rownames(xii$meta$cifti$labels[[1]])[4]
xii <- apply_xifti(xii, 1, function(x){ifelse(x==4, 1, 0)})
xii_d3 <- edit_mask_surf(
  xii$data$cortex_left[,1], mwall=rep(TRUE, 32492), 
  surf=load_surf("left"),
  do="dilate", depth=3
)
xii <- newdata_xifti(xii, cbind(as.matrix(xii), xii_d3))
xii <- convert_xifti(xii, "dlabel", colors="blue") # 17networks_LH_DefaultA_PCC_1

fnames <- paste0(c(tempfile(), tempfile(), tempfile()), ".png")
plot(xii, title=parc_name, idx=1, fname=fnames[1], legend_fname=fnames[3])
plot(xii, title="Dilated mask", idx=2, fname=fnames[2], legend_fname=fnames[3])
plt <- view_comp(fnames[seq(2)], fname=fnames[3], nrow=1)
knitr::include_graphics(plt, dpi=180)

500 evenly-sampled vertices

even_vert_samp obtains a subset of vertices sampled evenly across the mesh. It works by downsampling the mesh and then identifying the original vertices closest to the vertices on the downsampled mesh.

# Load a left-hemisphere surface.
surfL <- load_surf()

# Get 500 vertices evenly sampled on the mesh.
v <- even_vert_samp(surfL, 500)

# Get a mask of the 32k vertices indicating which are in the sample.
v2 <- seq(nrow(surfL$vertices)) %in% v

# Make a xifti having "1" for in-mask vertices, "0" elsewhere.
xii <- as.xifti(as.matrix(v2*1), surfL=surfL)

plot(convert_xifti(xii, "dlabel"))

500 evenly-sampled vertices

Loading included parcellations

ciftiTools includes several commonly-used parcellations which can be loaded in as a single-column dlabel "xifti" with load_parc. Refer to ?load_parc for their corresponding references. Any other parcellations can be read in with read_cifti.

parc <- load_parc("Schaefer_400")
# parc <- read_cifti("my_parcellation.dlabel.nii")

plot(parc)

## Warning in view_xifti_surface(xifti, color_mode = args$color_mode, zlim = args$zlim, : Too many labels (> 200) for qualitative color legend. Not rendering it.

Schaefer 400 parcellation

parc_add_subcortex will add new label levels and “parcels” data for each subcortical structure in the MNI template.

# Before `parc_add_subcortex`
max(parc) # 400 parcels

## [1] 400

summary(parc)

## ====CIFTI METADATA===================
## Intent:           3007 (dlabel)
## - names           "parcels"
## Measurements:     1 column
## 
## ====BRAIN STRUCTURES=================
## - left cortex     32492 data vertices
## 
## - right cortex    32492 data vertices

parc <- parc_add_subcortex(parc)
max(parc) # 400 cortex parcels + 19 subcortex parcels

## [1] 421

summary(parc) # After `parc_add_subcortex`

## ====CIFTI METADATA===================
## Intent:           3007 (dlabel)
## - names           "parcels"
## Measurements:     1 column
## 
## ====BRAIN STRUCTURES=================
## - left cortex     32492 data vertices
## 
## - right cortex    32492 data vertices
## 
## - subcortex       31870 data voxels
##                   subcortical structures and number of voxels in each:
##                     Cortex-L (0), Cortex-R (0),
##                     Accumbens-L (135), Accumbens-R (140),
##                     Amygdala-L (315), Amygdala-R (332),
##                     Brain Stem (3472),
##                     Caudate-L (728), Caudate-R (755),
##                     Cerebellum-L (8709), Cerebellum-R (9144),
##                     Diencephalon-L (706), Diencephalon-R (712),
##                     Hippocampus-L (764), Hippocampus-R (795),
##                     Pallidum-L (297), Pallidum-R (260),
##                     Putamen-L (1060), Putamen-R (1010),
##                     Thalamus-L (1288), Thalamus-R (1248),
##                     Other (0).

Applying a parcellation to data

apply_parc will apply a function to one "xifti" across each parcel specified by another "xifti". A common example is computing the average timeseries of each parcel:

parc <- load_parc("Yeo_7")
parc_tab <- parc$meta$cifti$labels[[1]]
rbind(head(parc_tab), tail(parc_tab)) # Keys 0-51 (52 total)

##                               Key      Red     Green      Blue Alpha
## ???                             0 1.000000 1.0000000 1.0000000     0
## 7Networks_LH_Vis                1 0.470588 0.0705882 0.5254900     1
## 7Networks_LH_SomMot             2 0.274510 0.5098040 0.7058820     1
## 7Networks_LH_DorsAttn_Post      3 0.000000 0.4627450 0.0549020     1
## 7Networks_LH_DorsAttn_FEF       4 0.000000 0.4627450 0.0588235     1
## 7Networks_LH_DorsAttn_PrCv      5 0.000000 0.4666670 0.0549020     1
## 7Networks_RH_Cont_PFCmp        46 0.901961 0.5803920 0.1411760     1
## 7Networks_RH_Default_Par       47 0.803922 0.2431370 0.3019610     1
## 7Networks_RH_Default_Temp      48 0.803922 0.2392160 0.3058820     1
## 7Networks_RH_Default_PFCv      49 0.803922 0.2392160 0.3019610     1
## 7Networks_RH_Default_PFCdPFCm  50 0.800000 0.2431370 0.3058820     1
## 7Networks_RH_Default_pCunPCC   51 0.800000 0.2431370 0.3019610     1

dat <- read_cifti(ciftiTools.files()$cifti["dtseries"])
dim(dat)

## [1] 60951     2

pmean <- apply_parc(
  dat, parc, 
  FUN=colMeans, # get average timeseries
  na.rm=TRUE # medial wall vals in `dat` will be NA
)
dim(pmean) # 52 parcels x 2 data columns

## [1] 52  2

For real timeseries data, cor(t(pmean)) would give the parcels x parcels functional connectivity matrix.

apply_parc also has an argument which converts the result to a "xifti" object. Each locations’ value will be the value of its corresponding parcel.

# Obtain the min/max, quartile, and mean value
#   of the second data column, within each parcel.
psumm <- apply_parc(
  select_xifti(dat, idx=2),
  parc, mwall_value=0, 
  FUN=summary,
  return_as="xifti"
)

# Plot the median.
plot(psumm, idx=3, title=psumm$meta$cifti$names[3], colors="magma")

## `zlim` not provided: using color range 0 - 3.2 (data limits: 0 - 3.41).

Median parcel value